#include <stdio.h>
#include "sep.h"
#include "septools.h"

typedef struct {
  char type;
  int nmols;
  int nmols_type;
} type_info;


double get_types(type_info *ptr, const char type, const char file[]);
void read_entry_from_file(double **omega, int i_entry, 
			  const char *file, int dir, type_info tinfo);

    

int main(int argc, char **argv){
  type_info tinfo;
  int flag = 1;
  int file_counter=0; 
  int nsample = 0;

  if ( argc != 4 ){
    fprintf(stderr, "Arguments: <type> <numb. of blocks per sample> <dir>\n"); 
    exit(EXIT_FAILURE);
  }

  char type = argv[1][0];
  int num_blocks_per_sample = atoi(argv[2]);
  int dir = atoi(argv[3]);

  double dt = get_types(&tinfo, type, "00001_mol.conf");
  
  double **omega = sep_matrix(num_blocks_per_sample, tinfo.nmols_type);
  double *acf = sep_vector(num_blocks_per_sample);

  while (1) {

    printf("\r");
    printf("Processing sample %d  ", nsample);
    fflush(stdout);

    for ( int n=0;  n<num_blocks_per_sample; n++ ){
      char file[256];
      sprintf(file, "%05d_mol.conf", file_counter);
      file_counter++;
      if ( file_exists(file) == 1 ){
	read_entry_from_file(omega, n, file, dir, tinfo);
	flag = 1;
      }
      else {
	flag = 0;
	break;
      }
    }
    
    if ( flag == 1 ){
      for ( int i=0; i<tinfo.nmols_type; i++ )
	for ( int k=0; k<num_blocks_per_sample; k++ )
	  for ( int kk=0; kk<num_blocks_per_sample-k; kk++ )
	    acf[k] += omega[kk][i]*omega[k+kk][i];   
      
      nsample++;
    }
    else if ( flag == 0 ) break;

  }
  printf("\n");

  FILE *fout = fopen("single_mavacf.dat", "w");
  if ( fout == NULL )
    sep_error("Couldn't open file - BAILING OUT");

  for ( int k=0; k< num_blocks_per_sample; k++ ){
    double t   = dt*(k+0.5);
    double fac = 1.0/((num_blocks_per_sample-k)*tinfo.nmols_type*nsample);
    fprintf(fout, "%f %f\n", t, acf[k]*fac);
  }

  fclose(fout);

  sep_free_matrix(omega, num_blocks_per_sample);
  free(acf);

  return 0;
}


double get_types(type_info *ptr, const char type, const char file[]){
  double tnow, lx, ly, lz;
  char line[2048];

  ptr->nmols = 0; ptr->nmols_type = 0; ptr->type = type;

  FILE *fin = fopen(file, "r");
  if ( fin == NULL ){
    sep_error("Couldn't read input file- BAILING OUT");
  }

  fscanf(fin, "%lf %lf %lf %lf\n", &tnow, &lx, &ly, &lz);
  
  while ( fgets(line, 2048, fin) != NULL ) {
    (ptr->nmols)++;
    if ( line[0] == type ) (ptr->nmols_type)++;
  }
  
  fclose(fin);

  printf("Found %d molecules of type %c (total number of molecules where %d)\n",
	 ptr->nmols_type, type, ptr->nmols);
  fflush(stdout);

  return tnow;

}

void read_entry_from_file(double **omega, int i_entry, 
			  const char *file, int dir, type_info tinfo){
  double tnow, lx, ly, lz;
  char type;
  int nuau;
  double m, x[3], v[3], s[3], w[3], inertia[6]; // xx,yy,zz,xy,xz,yz component
  
  FILE *fin = fopen(file, "r");
  if ( fin == NULL )
    sep_error("Couldn't read input file- BAILING OUT");
  
  fscanf(fin, "%lf %lf %lf %lf\n", &tnow, &lx, &ly, &lz);
  
  int counter = 0;
  for ( int i=0; i<tinfo.nmols; i++){
    int nread = fscanf(fin, "%c %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf\n", 
		       &type, &nuau, &m, 
		       &x[0], &x[1], &x[2], 
		       &v[0], &v[1], &v[2],
		       &s[0], &s[1], &s[2],
		       &w[0], &w[1], &w[2],
		       &inertia[0], &inertia[1], &inertia[2],
		       &inertia[3], &inertia[4], &inertia[5]);
    if ( nread != 21 )
      sep_error("Read failure - BAILING OUT");

    // Only this needs to be changed for other quantities
    if ( type == tinfo.type ){
      omega[i_entry][counter] = w[dir];
      counter++;
    }

  }
  
  fclose(fin);
} 
 

